home *** CD-ROM | disk | FTP | other *** search
/ Atari Mega Archive 1 / Atari Mega Archive - Volume 1.iso / apps / math / ols.zoo / olse.c < prev    next >
C/C++ Source or Header  |  1993-04-16  |  4KB  |  160 lines
  1. /*
  2.     Engine for doing OLS.
  3. */
  4.  
  5. #include <stdio.h>
  6. #include <math.h>
  7. #include <stddef.h>
  8. #include "utils.h"        /* for allocation/deallocation stuff. */
  9. #include "matrix.h"        /* for matrix inversion, dot products, etc. */
  10. #include "linstats.h"        /* for error checks. */
  11.  
  12. int 
  13.   olsengine (int noinference, float **data, int K, int ObsNum,
  14.          double *beta, double *S,
  15.          double *sigma2, double *R2, double *F,
  16.          float *Yhat)
  17.  
  18.      /*
  19.    Blackbox which does OLS.  Takes data matrix, K and ObsNum.
  20.    Produces vector of coefficients beta, vector of standard
  21.    errors S and scalar Rsquared.
  22.  
  23.    You have to allocate everything (beta, S and Yhat) before
  24.    calling ols.  Notice R2 and sigma2 are passed by reference.
  25.  
  26.    This function merely computes things; No I/O takes place.
  27. */
  28.  
  29. {
  30.   double **XpX;            /* X'X matrix. */
  31.   double **Xpy;            /* X'y matrix */
  32.   double **V;            /* (X'X) inverse matrix. */
  33.   double **tmpbeta;        /* temporary representation for beta. */
  34.   int i;
  35.  
  36.   if ((XpX = dmatrix (1, K, 1, K)) == NULL)
  37.     {
  38.       printf ("No room for X'X matrix.\n");
  39.       return 1;
  40.     }
  41.   if ((V = dmatrix (1, K, 1, K)) == NULL)
  42.     {
  43.       printf ("No room for X'X inverse matrix.\n");
  44.       return 1;
  45.     }
  46.   if ((Xpy = dmatrix (1, K, 1, 1)) == NULL)
  47.     {
  48.       printf ("No room for X'y vector.\n");
  49.       return 1;
  50.     }
  51.   if ((tmpbeta = dmatrix (1, K, 1, 1)) == NULL)
  52.     {
  53.       printf ("No room for beta vector.\n");
  54.       return 1;
  55.     }
  56.  
  57.   MakeXpX (data, XpX, K, ObsNum);
  58.   MakeXpy (data, Xpy, K, ObsNum);
  59.   if (1 == d_Inverse (XpX, V, K))
  60.     return 1;
  61.   d_MMult (V, Xpy, tmpbeta, K, K, 1);
  62.   /* multiply V (KxK) by Xpy (Kx1) giving tmpbeta (Kx1) */
  63.  
  64.   /* Now move into vector representation of beta: */
  65.   for (i = 0; i < K; i++)
  66.     beta[i] = tmpbeta[i + 1][1];
  67.   if (noinference)
  68.     return 0;
  69.   /* don't waste time doing inference when he says
  70.         he doesn't want it. */
  71.  
  72.   /* Generate predictions: */
  73.   Predict (data, beta, Yhat, K, ObsNum);    /* produces Yhat */
  74.   ANOVA (data, Yhat, K, ObsNum, sigma2, R2, F);
  75.   for (i = 0; i < K; i++)
  76.     S[i] = sqrt ((*sigma2) * V[i + 1][i + 1]);
  77.  
  78.   free_dmatrix (tmpbeta, 1, K, 1, 1);
  79.   free_dmatrix (Xpy, 1, K, 1, 1);
  80.   free_dmatrix (XpX, 1, K, 1, K);
  81.   free_dmatrix (V, 1, K, 1, K);
  82.  
  83.   return 0;
  84. }
  85.  
  86. void 
  87. Predict (float **data, double *beta, float *Yhat,
  88.      int K, int ObsNum)
  89. {
  90.   int i;
  91.  
  92.   for (i = 0; i < ObsNum; i++)
  93.     Yhat[i] = (float) df_dot (beta, data[i], K);
  94. }
  95.  
  96. void 
  97. MakeXpX (float **data, double **XpX, int K, int ObsNum)
  98. {
  99.   int vr, vc, i;
  100.   double sum;
  101.  
  102.   for (vr = 0; vr < K; vr++)
  103.     {
  104.       for (vc = 0; vc <= vr; vc++)
  105.     {
  106.       sum = 0.0;
  107.       for (i = 0; i < ObsNum; i++)
  108.         sum += (double) (data[i][vc] * data[i][vr]);
  109.       XpX[vr + 1][vc + 1] = sum;
  110.       XpX[vc + 1][vr + 1] = sum;
  111.     }
  112.     }
  113. }
  114.  
  115. void 
  116. MakeXpy (float **data, double **Xpy, int K, int ObsNum)
  117. {
  118.   int v, i;
  119.   double sum;
  120.  
  121.   for (v = 0; v < K; v++)
  122.     {
  123.       sum = 0.0;
  124.       for (i = 0; i < ObsNum; i++)
  125.     sum += (double) (data[i][v] * data[i][K]);
  126.       Xpy[v + 1][1] = sum;
  127.     }
  128. }
  129.  
  130. void 
  131. ANOVA (float **data, float *Yhat, int K, int ObsNum,
  132.        double *sigma2, double *R2, double *F)
  133.      /* Basically does a analysis of variance and reports results.*/
  134. {
  135.   double sumY, mean, dev;
  136.   int i;
  137.   double expl_var, total_var, sse;
  138.   /* explained variation, total variation and SSE */
  139.  
  140.   for (sumY = 0.0, i = 0; i < ObsNum; i++)
  141.     sumY += (double) data[i][K];
  142.   mean = sumY / ObsNum;
  143.  
  144.   expl_var = 0.0;
  145.   total_var = 0.0;
  146.   sse = 0.0;
  147.   for (i = 0; i < ObsNum; i++)
  148.     {
  149.       dev = (double) (mean - data[i][K]);
  150.       total_var += (dev * dev);
  151.       dev = (double) (mean - Yhat[i]);
  152.       expl_var += (dev * dev);
  153.     }
  154.   sse = total_var - expl_var;
  155.  
  156.   *sigma2 = sse / (ObsNum - K);
  157.   *R2 = expl_var / total_var;
  158.   *F = expl_var / ((K - 1) * (*sigma2));
  159. }
  160.